home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Interactive Reference Guide / C-C++ Interactive Reference Guide.iso / c_ref / csource3 / 188_01 / trans.c < prev    next >
Encoding:
C/C++ Source or Header  |  1987-10-01  |  6.9 KB  |  217 lines
  1. /*
  2. HEADER:        ;
  3. TITLE:        C elementary transcendentals;
  4. VERSION:    1.0;
  5.  
  6. DESCRIPTION:       Source    code    for   all   standard    C 
  7.  
  8.                transcendentals
  9.  
  10.         Employs  ldexp()  and  frexp()  functions;  if 
  11.  
  12.                suitable  versions  of these are  not  provided 
  13.  
  14.                by  a given compiler,  the versions provided in 
  15.  
  16.                source  code  will require  adaptation  to  the 
  17.  
  18.                double float formats of the compiler.
  19.  
  20. KEYWORDS:    Transcendentals, library;
  21. SYSTEM:        CP/M-80, V3.1;
  22. FILENAME:    TRANS.C;
  23. WARNINGS:      frexp()   and  ldexp()  are   implementation 
  24.  
  25.                dependent.   The  compiler  employed  does  not 
  26.  
  27.                support    minus   (-)   unary   operators   in 
  28.  
  29.                initializer  lists,  which are required by  the 
  30.  
  31.                code.
  32. AUTHORS:    Tim Prince;
  33. COMPILERS:    MIX C, v2.0.1;
  34. */
  35. /* Copyright (C) 1986, Chandler Software,
  36. ** 4456 W. Maple Rd., Birmingham MI 48010-1923
  37. ** (313) 855-4587
  38. ** C adaptation of 1984 FORTRAN code
  39. ** prepared for non-profit distribution by C Users' Group */
  40. /*$NESTCMNT*/
  41. # include "a:stdio"
  42. /* dblfmt describes floating point format
  43. ** machine dependencies in case frexp() and ldexp() are not
  44. ** supplied */
  45. union dblfmt{
  46.   double flt;
  47.   struct{
  48.     char mant[7]; /* full reverse byte order
  49. ** IEEE int sign:1;unsigned expn:11;unsigned mant:4
  50. ** with DEC reverse byte pairing, can't cross 16-bit
  51. ** boundaries */
  52.     char expn;
  53.   }fld;    /* char same as unsigned :8 in this C */
  54. }
  55. #define ROUND(x) (int)((x>=0)-.5+x) /* Pascal ROUND */
  56. #define ODD(i) ((i)&1) /* also like Pascal */
  57. /* teach the preprocessor some algebra */
  58. #define PI 3.1415926535897932385
  59. #define HPI 1.5707963267948966192
  60. #define LN2 .69314718055994530942
  61. #define L2E 1.442695040888963407
  62. #define cos(x) sin(HPI-(x)) /* double sin() */
  63. #define tan(x) sin(x)/cos(x) /* double sin() */
  64. /* for future use
  65. #define fabs(x) (x>=0?x:-(x))
  66. /* double sin() */
  67. #define atan2(x,y) (y>0?atan((x)/(y)):y==0?(x>=0?HPI:-HPI):(x>=0?PI+atan((x)/(y)):atan((x)/(y))-PI))
  68. #define asin(x) atan2(x,sqrt((1-(x))*(1+x))) /* double sin(),sqrt() */
  69. #define acos(x) atan2(sqrt((1-(x))*(1+x)),x) /* double sin(),sqrt() */
  70. #define log(x) log2(x)*LN2 /* double log2() */
  71. #define exp(x) exp2((x)*L2E) /* double exp2() */
  72. #define cosh(x) sqrt(sinh(x)*sinh(x)+1) /* double sinh(),sqrt() */
  73. #define tanh(x) (x<-20?-1.:(x>=20?1.:sinh(x)/cosh(x)))
  74. */
  75. #define pow(x,y) exp2(log2(x)*(y))
  76. main(){ /* test accuracy of tan,sin,cos,atan,exp2,log2 */
  77.   double sin(),atan(),log2(),exp2();
  78.   double x;
  79.   for(x=1e-36;x<1e33;x*=1e4)
  80.     printf("x:%23.16e tan(atan):%23.16e \n\tpow:%23.16e\n",
  81.       x,tan(atan(x)),pow(x,1.));
  82. }
  83. double frexp(x,scale)
  84.   int *scale;
  85.   double x;
  86. {
  87.   union dblfmt x1;
  88.   #define BIAS 128 /* exponent field of .9 (IEEE 1023) */
  89.   x1.flt=x;
  90.   *scale=x1.fld.expn-BIAS; /* scale by .5^(*scale) */
  91.   x1.fld.expn=BIAS;  /* .5 <= result < 1 */
  92.   return(x1.flt);
  93. }
  94. double ldexp(x,scale)
  95.   double x;
  96.   int scale;
  97. {
  98.   union dblfmt x1;
  99.   x1.flt=x;
  100. /* scale by 2^scale */
  101.   x1.fld.expn+=scale;
  102.   return(x1.flt);
  103. }
  104. double sin(x)
  105.   double x;
  106. {
  107.   #define NTS 8
  108. /* negative signs are ignored by some compilers! */
  109.   static double table[NTS]={
  110.     -.7370812e-12,.160478576e-9,-.2505187133e-7,
  111.     .275573164289e-5,-.00019841269823293,
  112.     .008333333333276252,-.1666666666666597,.9999999999999999};
  113.   double poly(),pm,floor();
  114. /* try to take care of x>maxint unless much too slow
  115. ** employ periodicity sin(x+n*PI)=(-1)^n*sin(x) */
  116.   x-=PI*(pm=floor(x/PI+.5));
  117. /* now |x| < HPI; use series and fix sign */
  118.   return(poly(x*x,NTS,table)*(ODD((int)pm)?-x:x));
  119. }
  120. /* use fast polynomial evaluation (possibly non-portable) */
  121. double poly(x,n,table)
  122.   double x;
  123.   int n;
  124.   double table[];
  125. {
  126.   register double sum;
  127.   register int i;
  128.   sum=table[0];
  129. /* unroll loop by pairing terms to save overhead */
  130.   for(i=2;i<n;i+=2)
  131.     sum=(sum*x+table[i-1])*x+table[i];
  132.   return(i==n?sum*x+table[i-1]:sum);
  133. }
  134. double atan(x)
  135.   double x;
  136. {
  137.   #define TNPI6 .57735026918962576451
  138.   #define NTA 9
  139.   static double table[NTA]={
  140.     .0443911883527,-.06483241134718,.07679374240257,
  141.     -.090903714847573,.111110979102203,-.1428571410307564,
  142.     .1999999999873222,-.33333333333329944,
  143.     .99999999999999999};
  144.   char invert,reduce,neg;
  145.   double poly();
  146. /* atan(-x)=-atan(x) */
  147.   if(neg=(x<0))x=-x;
  148. /* tan(HPI-x)=1/tan(x) */
  149.   if(invert=(x>=1))x=1/x;
  150. /* tan(a-b)={tan(a)-tan(b)}/{tan(a)*tan(b)+1}; b=PI/6 */
  151.   if(reduce=(x>=.269))x=(x-TNPI6)/(x*TNPI6+1);
  152.   x*=poly(x*x,NTA,table);
  153.   if(reduce)x+=.52359877559829887308;
  154.   if(invert)x=HPI-x;
  155.   return(neg?-x:x);
  156. }
  157. double log2(x)
  158.   double x;
  159. {
  160.   #define NTL 7
  161.   static double table[NTL]={
  162.     .24291838162,.2614440423316,.3206163946133,.41219840197457,
  163.     .57707801724629,.961796693924339,2.885390081777927};
  164.   int scale,incsc;
  165.   double poly(),ldexp(),frexp();
  166. /* split x -> 2^scale*(reduced x near 1) */
  167.   incsc=frexp(x,&scale)<.7071;
  168.   x=ldexp(x,-(scale-=incsc));
  169.   x=(x-1)/(x+1);
  170.   return(poly(x*x,NTL,table)*x+scale);
  171. }
  172. double exp2(x)
  173.   double x;
  174. {
  175.   #define INFINITY 1.7e38/*IEEE 0x7ff0000000000000*/
  176.   #define MAXEXP 127 /* IEEE 1023 */
  177.   #define MINLG2 -128
  178.   double evenp,oddp,xsq,ldexp();
  179.   int scale;
  180.   if(x>=MAXEXP)return(INFINITY);
  181.   if(x<MINLG2)return(0);
  182. /* 2^x=2^ROUND(x)*2^(x-ROUND(x)) */
  183.   x-=scale=ROUND(x);
  184. /* approximation is ratio of polynomials differing only in
  185. ** sign of odd terms */
  186.   xsq=x*x;
  187. /* continued fraction approx for e^x
  188.   (x*(15120+xsq*(420+xsq))+30240+xsq*(3360+xsq*30))/... */
  189. /* minimax fit for 2^x, 18 significant digits */
  190.   oddp=x*(65556.325877407443+xsq*(874.804294426022+xsq));
  191.   evenp=189155.572484473087+xsq*(10097.515376265486+
  192.     xsq*43.302654542155);
  193.   return(ldexp((evenp+oddp)/(evenp-oddp),scale));
  194. }
  195. /* for future use
  196. double sinh(x)
  197.   double x;
  198. {
  199.   #define NTSH 7
  200.   static double table[NTSH]={
  201.     1.632721e-10,2.50484370e-8,2.7557344083e-6,
  202.     1.984126975507e-4,.0083333333334753,.1666666666666579,
  203.     1.0000000000000001};
  204.   double poly(),exp2();
  205.   return(fabs(x)<1?poly(x*x,NTSH,table)*x:(exp(x)-1/exp(x))/2);
  206. }
  207. double sqrt(x) /* as if division hardware but no sqrt */
  208.   double x;
  209. { static float table[]={.59,.417,.295}; /* minimax 2 digits */
  210.   register double x1;
  211.   double frexp(),ldexp();
  212.   register int i;
  213.   register char iters=3; /* number of Newton iterations */
  214.   register float x0; /* early stages need 2.5 digits precision */
  215.   int scale;
  216.   if(x<=0)return(x);
  217.   x0=frexp(x,&scale); /* sqrt(2^x*y)=2^(x/2)sqrt(y) */
  218.   i=ODD(scale); /* separate fits for y<.5 & y>.5 */
  219.   x1=ldexp(x0*table[i]+table[i+1],(scale+1)>>1); /*crude sqrt*/
  220.   do /* enough Newton iterations for full accuracy */
  221.     x1=ldexp(x/x1+x1,-1); /* ldexp might be faster than *.5 */
  222.   while(--iters);
  223.   return(x1);
  224. }
  225. */
  226.